Progetto Bioinformatica Rocco Saviano
Rocco Saviano
2024-05-31
Analisi dell’espressione differenziale in diverse condizioni
Nella seguente analisi, l’obiettivo è studiare l’espressione differenziale di specifici geni in diverse condizioni, in particolare nel contesto del cancro al seno. Utilizziamo esemplari di topo, che condividono gran parte del genoma con l’essere umano, per studiare l’espressione differenziale in tre condizioni: silenziamento di BRCA1 (gene che produce una proteina coinvolta nella riparazione del DNA; la sua disfunzione è associata a un aumentato rischio di sviluppare certi tipi di cancro, come il cancro al seno), silenziamento di BRCA2 (gene che produce una proteina fondamentale per la riparazione del DNA e la stabilità genomica. Le proteine codificate da questi geni aiutano a riparare i danni al DNA che possono portare a mutazioni e sviluppo di tumori), e la condizione “Wild Type” (senza silenziamento).
Analisi Silenziamento BRCA2 vs WT senza silenziamento
In questo primo confronto si intende confrontare l’espressione differenziale tra i campioni tumorali che presentano un silenziamento del gene BRCA2 e i campioni tumorali che non presentano nessun tipo di silenziamento
Step 1: Import dei dati
Sfruttiamo la libreria GEOquery per caricare i dati relativi a questo specifico esperimento e successivamente consideriamo solo la matrice dei conteggi grezza di interesse da cui estraiamo i conteggi
library(GEOquery)
gse <- getGEO("GSE137818")[[1]]
#getGEOSuppFiles("GSE137818")
raw_data <- read.csv("GSE137818/GSE137818_Mouse_UNTREATED_bulkRNA_BRCA2KO_vs_WT_rawcounts.csv.gz", stringsAsFactors = FALSE)
head(raw_data[,1:5])## X KO_d14_IgG_1 KO_d14_IgG_2 KO_d14_IgG_3 KO_d14_IgG_4
## 1 100009600 23 29 57 28
## 2 100009614 0 0 0 0
## 3 100009664 1 1 0 0
## 4 100017 2330 2586 2683 2705
## 5 100019 3428 3522 2736 3616
## 6 100033459 217 265 166 101conteggi <- raw_data[,-1]Step 2: Modifica e processing dei dati
In questa fase andiamo a modificare i nomi dei campioni rispetto al formato orginale per migliorarne l’interpretazione
old_names <- strsplit(names(raw_data), split="_")
new_names <- vector("list", length(old_names))
for (i in 2:length(old_names)) {
if (old_names[[i]][1] == "KO") {
old_names[[i]][1] <- "BRCA2"
}
if (length(old_names[[i]]) > 1) {
new_names[[i]] <- paste(old_names[[i]][1], old_names[[i]][length(old_names[[i]])], sep = "-")
} else {
new_names[[i]] <- old_names[[i]][1]
}
}
new_names <- unlist(new_names)
names(conteggi) <- new_names #abbiamo cambiato i nomi rendendoli un po' più leggibili
head(conteggi)## BRCA2-1 BRCA2-2 BRCA2-3 BRCA2-4 BRCA2-5 WT-1 WT-2 WT-3 WT-4 WT-5
## 1 23 29 57 28 28 16 24 23 14 10
## 2 0 0 0 0 0 0 0 0 0 0
## 3 1 1 0 0 0 0 0 0 0 0
## 4 2330 2586 2683 2705 3069 1659 1978 1973 1437 1591
## 5 3428 3522 2736 3616 3941 2353 3357 2990 2366 3362
## 6 217 265 166 101 119 36 153 77 59 52Andiamo a costruire l’oggetto della classe SummarizedExpreriment, ma prima andiamo a costruire i rowData formati dalla coppia ENTREZID-GENE SYMBOL per ogni gene, i colData formati semplicemente dai nomi dei campioni e dalla condizione Silenziato/non silenziato
library(SummarizedExperiment)
library(clusterProfiler)
id_geni <- raw_data[,1]
rownames(conteggi) <- id_geni
organism <- 'org.Mm.eg.db'
gene_symb <- bitr(id_geni, fromType="ENTREZID", toType="SYMBOL",OrgDb = organism)
row_data <- data.frame(ENTREZID=id_geni, ordine=seq_along(id_geni))
row_data <- merge(row_data ,gene_symb, by = "ENTREZID", all.x=TRUE)Abbiamo poi costruito i rowRanges cioè le posizioni dei geni sul cromosoma di riferimento e il cromosoma di riferimento, per fare ciò abbiamo sfruttato il pacchetto biomaRt che ci permette di collegarci ad un database di Ensembl in particolare a quello del topo, successivamente abbiamo effettuato il mapping grazie al pacchetto GenomicRanges. Infine è stato possibile costruire il summarized Experiment contente tutte le informazioni neccessarie per l’analisi, andiamo poi a salvare il SummarizedExperiment in modo da poterlo richiamare facilmente nelle prossime occasioni
library(dplyr)
library(biomaRt)
library(GenomicRanges)
fonte <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")
posizioni_geni <- getBM(filters = "entrezgene_id", attributes = c("entrezgene_id", "chromosome_name", "start_position", "end_position"), values = id_geni, mart = fonte)
names(posizioni_geni)[1] <- "ENTREZID"
#nel caso di duplicati teniamo in considerazione solo la prima occorrenza
posizioni_geni <- posizioni_geni %>% distinct(ENTREZID, .keep_all = TRUE)
row_data <- merge(row_data, posizioni_geni, by = "ENTREZID", all.x=TRUE)
row_data <- row_data[order(row_data$ordine),]
row_data$ordine <- NULL
silenziamento <- lapply(new_names, function(name){
if (substr(name, 1, 4) == "BRCA") {
return("S")
} else {
return("NS")
}
})
silenziamento <- unlist(silenziamento)
colData <- cbind(new_names, silenziamento)
se <- SummarizedExperiment(
assays = conteggi,
rowData = row_data,
colData = colData,
)head(assay(se))## BRCA2-1 BRCA2-2 BRCA2-3 BRCA2-4 BRCA2-5 WT-1 WT-2 WT-3 WT-4 WT-5
## 100009600 23 29 57 28 28 16 24 23 14 10
## 100009614 0 0 0 0 0 0 0 0 0 0
## 100009664 1 1 0 0 0 0 0 0 0 0
## 100017 2330 2586 2683 2705 3069 1659 1978 1973 1437 1591
## 100019 3428 3522 2736 3616 3941 2353 3357 2990 2366 3362
## 100033459 217 265 166 101 119 36 153 77 59 52colData(se)## DataFrame with 10 rows and 2 columns
## new_names silenziamento
## <character> <character>
## BRCA2-1 BRCA2-1 S
## BRCA2-2 BRCA2-2 S
## BRCA2-3 BRCA2-3 S
## BRCA2-4 BRCA2-4 S
## BRCA2-5 BRCA2-5 S
## WT-1 WT-1 NS
## WT-2 WT-2 NS
## WT-3 WT-3 NS
## WT-4 WT-4 NS
## WT-5 WT-5 NSsaveRDS(se, file = "GSE137818_SummarizedExperiment.rds")
rowData(se)## DataFrame with 20637 rows and 5 columns
## ENTREZID SYMBOL chromosome_name start_position end_position
## <integer> <character> <character> <integer> <integer>
## 100009600 100009600 Zglp1 9 20973689 20978389
## 100009614 100009614 Gm10024 10 77547291 77547843
## 100009664 100009664 F630206G17Rik NA NA NA
## 100017 100017 Ldlrap1 4 134468865 134495335
## 100019 100019 Mdn1 4 32657119 32775217
## ... ... ... ... ... ...
## 99889 99889 Arfip1 3 84403400 84489932
## 99890 99890 Prmt6 3 110153425 110158314
## 99899 99899 Ifi44 3 151436559 151455597
## 99929 99929 Tiparp 3 65435831 65462939
## 99982 99982 Kdm1a 4 136277851 136330034sort(table(rowData(se)$chromosome_name))##
## GL456221.1 Y 18 16 19 12 X
## 1 2 521 636 660 683 703
## 14 13 15 10 17 8 3
## 711 777 805 954 988 1008 1010
## 6 9 1 5 4 2 7
## 1050 1101 1190 1208 1234 1556 1580
## 11
## 1612se <- readRDS("GSE137818_SummarizedExperiment.rds")Step 3: Analisi esplorativa dei dati e filtraggio
Possiamo fare un primo boxplot su scala logaritmica per capire se i dati che guardiamo sono comparabili
library(RColorBrewer)
pal<-brewer.pal(3, "Set2")
boxplot(log1p(assay(se)), col=pal[as.factor(se$silenziamento)], las=2)
e osserviamo che tutto sommato per tutti i campioni siamo sugli stessi ordini di grandezza, possiamo di seguito effettuare un filtraggio andando a scartare i geni che sono poco presenti in particolare quelli con somma di riga minore o uguale a 10
sef <- se[rowSums(assay(se))>10,]
assay(sef,"log") <- log1p(assay(sef))
boxplot(assay(sef,"log"), col=pal[as.factor(sef$silenziamento)], las=2)
A seguito del filtraggio possiamo osservare che rimuovendo i geni scarsamente espressi che potrebbero rappresentare un rumore nei nostri dati, abbiamo migliorato le distribuzioni dei dati perchè è come se avessimo delle distribuzioni più “regolari” e simmetriche rispetto al boxplot precedente
Andiamo poi ad fare un grafico RLE (Relative Log Expression) che ci restituisce come informazione il logaritmo del valore dei nostri geni normalizzati rispetto alla loro mediana su tutti i campioni (log(gi_c1/gi_cmedian)), possiamo osservare che già senza nessuna normalizzazione dovremmo avere una lieve sovraespressione nei campioni con silenziamento di BRCA2 a meno di differenze legate all’atto dell’esperimento, effetto che proveremo poi a mitigare senza appiattire troppo il contenuto informativo del segnale successivamente con la normalizzazione.
library(EDASeq)
plotRLE(as.matrix(assay(sef)), col=pal[as.factor(sef$silenziamento)], outline=FALSE, las=2)
Andando poi a costruire la PCA osserviamo la netta separazione delle due tipologie di campione, anche se i campioni WT risultano poco vicini in termini della PC2 e la varianza totale spiegata è di circa 50% rispetto al 70-80% che ci aspettiamo per una buona rappresentazione dei dati
plotPCA(as.matrix(assay(sef)), col=pal[as.factor(sef$silenziamento)])
Risulta quindi necessario a valle di queste ultime osservazioni cercare di migliorare la situazione andando a normalizzare i dati
Step 4: Esplorazione delle normalizzazioni
Calcoliamo diverse tipologie di normalizzazione ovvero:
-Upper Quartile: Allinea i dati rispetto al valore del 3 quartile
-Full Quartile: Usa tutti i quartili per calcolare un fattore di scala
-Trimmed Mean of M values: Calcola un fattore di scala sfruttando la media dei rapporti di espressione dei geni
-RLE: “median library is calculated from the geometric mean of all columns and the median ratio of each sample to the median library is taken as the scale factor.”
library(edgeR)
assayNames(sef)[1] <- "counts"
assay(sef, "upper") <- betweenLaneNormalization(as.matrix(assay(sef)), which="upper")
assay(sef, "full") <- betweenLaneNormalization(as.matrix(assay(sef)), which="full")
assay(sef, "tmm") <- cpm(assay(sef), lib.size = calcNormFactors(assay(sef), method = "TMM") * colSums(assay(sef)))
assay(sef, "rle_gm") <- cpm(assay(sef), lib.size = calcNormFactors(assay(sef), method = "RLE") * colSums(assay(sef)))for(a in c("counts","upper", "full", "tmm", "rle_gm")) {
plotRLE(as.matrix(assay(sef, a)), col=pal[as.factor(sef$silenziamento)], las=2, outline=FALSE, main = a)
}
for(a in c("counts","upper", "full", "tmm", "rle_gm")) {
plotPCA(as.matrix(assay(sef, a)), col=pal[as.factor(sef$silenziamento)], main = a)
}
Step5: Stima e validazione dei parametri di normalizzazione e dispersione
A seguito dello step precedente, la normalizzazione che spiega la maggior parte della varianza, permette una netta separazione dei campioni in due regioni dello spazio delle componenti principali senza far collassare i diversi campioni in un unico punto e allo stesso tempo permette nei boxplot di avere mediane allineate e box confrontabili è la “TMM”. Andiamo quindi a calcolare i fattori di normalizzazione e li salviamo nella variabile ‘dge’
dge <- calcNormFactors(sef, method = "TMM")Andiamo poi a costruire la model matrix, cioè come le covariate (cioè le variabili indipendenti che nel nostro caso è il Silenziamento) sono organizzate (se più di 1, non come nel nostro caso) per modellare la loro dipendenza dalla variabile dipendente che nel nostro caso è l’espressione genica dei vari geni
design <- model.matrix(~ silenziamento ,data=colData(sef))Andiamo poi a stimare la dispersione dei dati grazie alla funzione “estimateDisp” che sfruttando i fattori di normalizzazione e la matrice di design ci restituisce:
-la dispersione globale comune a tutti i geni
-relazione tra media e dispersione di ogni singolo gene (come la varianza cambia in funzione del livello medio di espressione genica)
-la dispersione specifica per ogni gene
Successivamente andiamo a plottare il grafico “Media-Varianza” in cui:
-Gli elementi grigi rappresentano la varianza grezza cioè la varianza osservata per ciascun gene
-Gli elementi azzurri rappresentano la varianza stimata per ciascun gene sfruttando i parametri di dispersione calcolati nello stesso blocco di codice
-La trend-line rappresenta la varianza attesa dato un certo livello medio di espressione genica, considerando una distribuzione binomiale negativa
dge <- estimateDisp(dge, design)
plotMeanVar(dge, show.raw.vars = TRUE, show.tagwise.vars = TRUE,
show.ave.raw.vars = FALSE)
Possiamo quindi dedurre dal grafico Media-Varianza che sia i dati grezzi che quelli stimati presentano una varianza maggiore del normale (rispetto alla linea nera), suggerendo la presenza di geni con variazioni significative di espressione tra i vari campioni che potrebbero quindi essere differenzialmente espressi.
A seguire, andiamo a plottare il coefficiente di variazione biologica (BCV) che plotta i conteggi rispetto ad un parametro che misura la dispersione rispetto alla variabilità biologica dei campioni Nel grafico possiamo osservare:
-I pallini neri che rappresentano il coefficiente di variazione biologica stimata per ogni campione
-il coefficiente di variazione biologica comune stimato su tutti i geni
-il coefficiente di variazione biologica stimata rispetto al livello medio di espressione di tutti i geni
plotBCV(dge)
Possiamo a valle di questi grafici fare le seguenti osservazioni: Osserviamo che i tagwise seguono l’andamento della Trend per la maggior parte dei geni, inoltre osserviamo che seppur nel grafico Mean-Variance precedente i geni non fossero perfettamente distribuiti intorno alla varianza attesa, il trend viene seguito in maniera abbastanza fedele, dunque queste informazioni ci danno sicurezza rispetto alla “validazione” dei parametri di stima calcolati con la funzione ‘estimateDisp’
Step 6: Analisi dell’espressione differenziale
Con i parametri di dispersione (per catturare la variabilità biologica) e i coefficienti di normalizzazione (per eliminare variazioni non biologiche) ottenuti e validati nello step precedente, andiamo a costruire lo specifico modello lineare generalizzato (GLM) specificando il contrasto per indicare che nel calcolo del log(FC) debba considerare Silenzati/Non Silenziati e cambiando la matrice di design aggiungendo uno 0 in modo da usare un “non-intercept model” come di consueto in questo tipo di analisi.
design0 <- model.matrix(~ 0 + silenziamento ,data=colData(sef))
cont <- makeContrasts(silenziamentoS - silenziamentoNS, levels=design0)
glm_s <- glmFit(dge, design0)Eseguiamo il test di espressione differenziale utilizzando il modello costruito
ris <- glmLRT(glm_s, contrast = cont)top <- topTags(ris, n=Inf)$table
diff_exp <- top[top$FDR<=0.05,]
up_reg <- diff_exp[(diff_exp$logFC)>0,]
down_reg <- diff_exp[(diff_exp$logFC)<0,]
table(top$FDR<=0.05)##
## FALSE TRUE
## 12872 4273Ottenendo 4273 geni differenzialmente espressi considerando come p-value ajusted (FDR) ≤ 0.05. che si distribuiscono sui diversi cromosomi nel seguente modo:
sort(table(diff_exp$chromosome_name))##
## GL456221.1 18 16 19 12 13 X
## 1 89 129 137 138 139 141
## 14 17 9 6 8 10 15
## 143 185 202 203 204 205 214
## 4 5 2 3 7 11 1
## 243 263 289 297 300 327 341Inoltre abbiamo suddiviso il nostro set in due partizioni, i geni up-regolati che rappresentano una maggiore presenza nel caso di silenziamento di BRCA2 e i geni down-regolati che rappresentano una minore presenza nel caso di silenziamento di BRCA2
dim(up_reg)## [1] 2276 10dim(down_reg)## [1] 1997 10Step 7: Analisi funzionale ed enrchiment
Ottenuti i geni che risultano espressi in maniera significativamente diversa tra i due cluster di campioni, possiamo andare a fare un’analisi funzionale per cercare di capire a che servono questi geni e quindi in che modo la loro diversa espressione si ripercuote poi in termini funzionali
Analisi Generale
La prima cosa è andare a guardare in generale, senza distinzione di up e down regolati cosa possiamo dedurre
Gene Ontology
Andiamo a fare un’analisi di arricchimento rispetto alla Gene Ontology, e successivamente andiamo a fare diversi plot in diverse condizioni considerando per il cnetplot una partizione relativa ai primi 15 elementi relativi alla gene ontology per avere una maggiore leggibilità
library(simplifyEnrichment)
library(clusterProfiler)
library(enrichplot)
gene_sym <- diff_exp$SYMBOL
go_sym <- enrichGO(gene_sym, OrgDb = organism,
keyType ="SYMBOL", ont = "ALL")
dotplot(go_sym, showCategory=10)
go_sym1<-pairwise_termsim(go_sym)
emapplot(go_sym1, showCategory=10)
gene_list <- diff_exp
fold_changes <- setNames(gene_list$logFC, gene_list$SYMBOL)
go_sym_sub <- subset_enrichResult(go_sym, 15)
cnetplot(go_sym, foldChange=fold_changes, circular = TRUE, colorEdge = TRUE)
cnetplot(go_sym_sub, foldChange=fold_changes, circular = TRUE, colorEdge = TRUE)
Usando invece i plot offerti da gprofiler2
library(gprofiler2)
geneOnt_ris <- gost(gene_sym, sources="GO",organism = "mmusculus", significant = TRUE, user_threshold = 0.05)
gostplot(geneOnt_ris, interactive = TRUE)Andando però a dividere le tre macro aree del vocabolario gene ontology osserviamo che
table(go_sym@result$ONTOLOGY)##
## BP CC MF
## 2607 186 312((go_sym@result)[go_sym@result$ONTOLOGY == "BP", ])[1:5,3]## [1] "regulation of innate immune response"
## [2] "positive regulation of innate immune response"
## [3] "immune response-regulating signaling pathway"
## [4] "positive regulation of response to biotic stimulus"
## [5] "response to virus"((go_sym@result)[go_sym@result$ONTOLOGY == "CC", ])[1:5,3]## [1] "collagen-containing extracellular matrix"
## [2] "cell-substrate junction"
## [3] "focal adhesion"
## [4] "cell leading edge"
## [5] "membrane raft"((go_sym@result)[go_sym@result$ONTOLOGY == "MF", ])[1:5,3]## [1] "GTPase regulator activity"
## [2] "nucleoside-triphosphatase regulator activity"
## [3] "GTP binding"
## [4] "guanyl nucleotide binding"
## [5] "guanyl ribonucleotide binding"KEGG Pathways - Overview
Sfruttiamo poi i KEGG Pathways per continuare la nostra analisi fuzionale usando sia ClusterProfiler che Gprofiler2
gene_entrez <- diff_exp$ENTREZID
Enrich.KEG <- enrichKEGG(gene_entrez, organism="mmu")
Enrich.KEG1<-pairwise_termsim(Enrich.KEG)
dotplot(Enrich.KEG , showCategory=10)
emapplot(Enrich.KEG1, , showCategory=10)
cnetplot(Enrich.KEG, foldChange=fold_changes, circular = TRUE, colorEdge = TRUE) 
Kegg_res <- gost(gene_entrez, sources="KEGG", organism = "mmusculus", significant = TRUE,user_threshold = 0.05)
gostplot(Kegg_res, interactive = TRUE)Abbiamo trovato quindi i seguenti risultati per Gprofiler2 e ClusterProfiler
Kegg_res$result$term_name## [1] "Toxoplasmosis"
## [2] "Focal adhesion"
## [3] "Th17 cell differentiation"
## [4] "T cell receptor signaling pathway"
## [5] "Th1 and Th2 cell differentiation"
## [6] "Fluid shear stress and atherosclerosis"
## [7] "ECM-receptor interaction"
## [8] "MAPK signaling pathway"
## [9] "Fc gamma R-mediated phagocytosis"
## [10] "Rheumatoid arthritis"
## [11] "Yersinia infection"
## [12] "Chemokine signaling pathway"
## [13] "Leishmaniasis"
## [14] "TNF signaling pathway"
## [15] "Coronavirus disease - COVID-19"
## [16] "Pertussis"
## [17] "Efferocytosis"
## [18] "Human papillomavirus infection"
## [19] "Measles"
## [20] "Proteoglycans in cancer"
## [21] "PI3K-Akt signaling pathway"
## [22] "Chagas disease"
## [23] "Influenza A"
## [24] "NOD-like receptor signaling pathway"
## [25] "Choline metabolism in cancer"
## [26] "PD-L1 expression and PD-1 checkpoint pathway in cancer"
## [27] "Lipid and atherosclerosis"
## [28] "Toll-like receptor signaling pathway"
## [29] "C-type lectin receptor signaling pathway"
## [30] "Tuberculosis"
## [31] "Osteoclast differentiation"
## [32] "Axon guidance"
## [33] "Metabolic pathways"
## [34] "Cytokine-cytokine receptor interaction"
## [35] "Leukocyte transendothelial migration"
## [36] "Sphingolipid signaling pathway"
## [37] "Fc epsilon RI signaling pathway"
## [38] "Lysosome"
## [39] "Virion - Herpesvirus"
## [40] "Rap1 signaling pathway"
## [41] "Inflammatory bowel disease"
## [42] "Arrhythmogenic right ventricular cardiomyopathy"
## [43] "Amoebiasis"
## [44] "Regulation of actin cytoskeleton"
## [45] "Ras signaling pathway"
## [46] "Endocrine resistance"
## [47] "Pathways in cancer"Enrich.KEG$Description## [1] "Toxoplasmosis - Mus musculus (house mouse)"
## [2] "Focal adhesion - Mus musculus (house mouse)"
## [3] "Th17 cell differentiation - Mus musculus (house mouse)"
## [4] "T cell receptor signaling pathway - Mus musculus (house mouse)"
## [5] "Th1 and Th2 cell differentiation - Mus musculus (house mouse)"
## [6] "MAPK signaling pathway - Mus musculus (house mouse)"
## [7] "Fluid shear stress and atherosclerosis - Mus musculus (house mouse)"
## [8] "ECM-receptor interaction - Mus musculus (house mouse)"
## [9] "Fc gamma R-mediated phagocytosis - Mus musculus (house mouse)"
## [10] "Rheumatoid arthritis - Mus musculus (house mouse)"
## [11] "Chemokine signaling pathway - Mus musculus (house mouse)"
## [12] "TNF signaling pathway - Mus musculus (house mouse)"
## [13] "Yersinia infection - Mus musculus (house mouse)"
## [14] "Leishmaniasis - Mus musculus (house mouse)"
## [15] "Measles - Mus musculus (house mouse)"
## [16] "Efferocytosis - Mus musculus (house mouse)"
## [17] "Human papillomavirus infection - Mus musculus (house mouse)"
## [18] "Pertussis - Mus musculus (house mouse)"
## [19] "Proteoglycans in cancer - Mus musculus (house mouse)"
## [20] "PI3K-Akt signaling pathway - Mus musculus (house mouse)"
## [21] "Chagas disease - Mus musculus (house mouse)"
## [22] "Influenza A - Mus musculus (house mouse)"
## [23] "Cytoskeleton in muscle cells - Mus musculus (house mouse)"
## [24] "Toll-like receptor signaling pathway - Mus musculus (house mouse)"
## [25] "Choline metabolism in cancer - Mus musculus (house mouse)"
## [26] "PD-L1 expression and PD-1 checkpoint pathway in cancer - Mus musculus (house mouse)"
## [27] "Lipid and atherosclerosis - Mus musculus (house mouse)"
## [28] "Tuberculosis - Mus musculus (house mouse)"
## [29] "NOD-like receptor signaling pathway - Mus musculus (house mouse)"
## [30] "Cytokine-cytokine receptor interaction - Mus musculus (house mouse)"
## [31] "C-type lectin receptor signaling pathway - Mus musculus (house mouse)"
## [32] "Axon guidance - Mus musculus (house mouse)"
## [33] "Osteoclast differentiation - Mus musculus (house mouse)"
## [34] "Sphingolipid signaling pathway - Mus musculus (house mouse)"
## [35] "Leukocyte transendothelial migration - Mus musculus (house mouse)"
## [36] "Lysosome - Mus musculus (house mouse)"
## [37] "Fc epsilon RI signaling pathway - Mus musculus (house mouse)"
## [38] "Virion - Herpesvirus - Mus musculus (house mouse)"
## [39] "Rap1 signaling pathway - Mus musculus (house mouse)"
## [40] "Inflammatory bowel disease - Mus musculus (house mouse)"
## [41] "Ras signaling pathway - Mus musculus (house mouse)"
## [42] "Regulation of actin cytoskeleton - Mus musculus (house mouse)"
## [43] "Amoebiasis - Mus musculus (house mouse)"
## [44] "Endocrine resistance - Mus musculus (house mouse)"
## [45] "Arrhythmogenic right ventricular cardiomyopathy - Mus musculus (house mouse)"
## [46] "Alcoholic liver disease - Mus musculus (house mouse)"
## [47] "Central carbon metabolism in cancer - Mus musculus (house mouse)"
## [48] "B cell receptor signaling pathway - Mus musculus (house mouse)"
## [49] "HIF-1 signaling pathway - Mus musculus (house mouse)"
## [50] "Epstein-Barr virus infection - Mus musculus (house mouse)"
## [51] "Neurotrophin signaling pathway - Mus musculus (house mouse)"
## [52] "Primary immunodeficiency - Mus musculus (house mouse)"
## [53] "Glycerophospholipid metabolism - Mus musculus (house mouse)"
## [54] "RIG-I-like receptor signaling pathway - Mus musculus (house mouse)"
## [55] "Viral protein interaction with cytokine and cytokine receptor - Mus musculus (house mouse)"
## [56] "ErbB signaling pathway - Mus musculus (house mouse)"
## [57] "mTOR signaling pathway - Mus musculus (house mouse)"
## [58] "Asthma - Mus musculus (house mouse)"
## [59] "Cytosolic DNA-sensing pathway - Mus musculus (house mouse)"
## [60] "Protein processing in endoplasmic reticulum - Mus musculus (house mouse)"
## [61] "Estrogen signaling pathway - Mus musculus (house mouse)"
## [62] "Hematopoietic cell lineage - Mus musculus (house mouse)"
## [63] "Relaxin signaling pathway - Mus musculus (house mouse)"
## [64] "Apoptosis - Mus musculus (house mouse)"
## [65] "Malaria - Mus musculus (house mouse)"
## [66] "Adipocytokine signaling pathway - Mus musculus (house mouse)"
## [67] "Thyroid hormone signaling pathway - Mus musculus (house mouse)"
## [68] "Prostate cancer - Mus musculus (house mouse)"
## [69] "Hypertrophic cardiomyopathy - Mus musculus (house mouse)"
## [70] "Cell adhesion molecules - Mus musculus (house mouse)"
## [71] "AGE-RAGE signaling pathway in diabetic complications - Mus musculus (house mouse)"
## [72] "Hepatitis C - Mus musculus (house mouse)"
## [73] "Pyruvate metabolism - Mus musculus (house mouse)"
## [74] "Sphingolipid metabolism - Mus musculus (house mouse)"
## [75] "Salmonella infection - Mus musculus (house mouse)"
## [76] "JAK-STAT signaling pathway - Mus musculus (house mouse)"
## [77] "IL-17 signaling pathway - Mus musculus (house mouse)"
## [78] "Hepatitis B - Mus musculus (house mouse)"
## [79] "Tight junction - Mus musculus (house mouse)"
## [80] "VEGF signaling pathway - Mus musculus (house mouse)"
## [81] "Endocytosis - Mus musculus (house mouse)"
## [82] "NF-kappa B signaling pathway - Mus musculus (house mouse)"
## [83] "Phospholipase D signaling pathway - Mus musculus (house mouse)"
## [84] "Dilated cardiomyopathy - Mus musculus (house mouse)"
## [85] "Legionellosis - Mus musculus (house mouse)"
## [86] "Platinum drug resistance - Mus musculus (house mouse)"
## [87] "Insulin resistance - Mus musculus (house mouse)"
## [88] "Kaposi sarcoma-associated herpesvirus infection - Mus musculus (house mouse)"
## [89] "Synaptic vesicle cycle - Mus musculus (house mouse)"
## [90] "Glycosaminoglycan biosynthesis - keratan sulfate - Mus musculus (house mouse)"
## [91] "Hepatocellular carcinoma - Mus musculus (house mouse)"
## [92] "Small cell lung cancer - Mus musculus (house mouse)"
## [93] "Coronavirus disease - COVID-19 - Mus musculus (house mouse)"
## [94] "African trypanosomiasis - Mus musculus (house mouse)"
## [95] "Antigen processing and presentation - Mus musculus (house mouse)"
## [96] "Biosynthesis of amino acids - Mus musculus (house mouse)"
## [97] "Human immunodeficiency virus 1 infection - Mus musculus (house mouse)"Che sono coerenti tra loro, possiamo inoltre notare che come visto a lezione è presente il pathway per il Covid19 ma non possiamo dire se questo silenziamento in termini di KEGG Pathways vada ad influire positivamente o nega tivamente per il cancro al seno, infatti come possiamo notare è presente “Pathways in cancer” ([47] di Kegg_res) ma non sappiamo se in termini positivi o negativi
KEGG Pathways - Geni UP REGOLATI
upgene_entrez <- up_reg$ENTREZID
Kegg_res <- gost(upgene_entrez, sources="KEGG", organism = "mmusculus", significant = TRUE,user_threshold = 0.05)
Kegg_res$result$term_name## [1] "Th1 and Th2 cell differentiation"
## [2] "Th17 cell differentiation"
## [3] "Toxoplasmosis"
## [4] "NOD-like receptor signaling pathway"
## [5] "Influenza A"
## [6] "Measles"
## [7] "Epstein-Barr virus infection"
## [8] "T cell receptor signaling pathway"
## [9] "Leishmaniasis"
## [10] "Coronavirus disease - COVID-19"
## [11] "Inflammatory bowel disease"
## [12] "Staphylococcus aureus infection"
## [13] "Herpes simplex virus 1 infection"
## [14] "Cell adhesion molecules"
## [15] "Pertussis"
## [16] "Cytokine-cytokine receptor interaction"
## [17] "Chemokine signaling pathway"
## [18] "Cytosolic DNA-sensing pathway"
## [19] "Toll-like receptor signaling pathway"
## [20] "Rheumatoid arthritis"
## [21] "Osteoclast differentiation"
## [22] "Asthma"
## [23] "Yersinia infection"
## [24] "Chagas disease"
## [25] "Primary immunodeficiency"
## [26] "Fc gamma R-mediated phagocytosis"
## [27] "Tuberculosis"
## [28] "Intestinal immune network for IgA production"
## [29] "RIG-I-like receptor signaling pathway"
## [30] "Hepatitis C"
## [31] "Leukocyte transendothelial migration"
## [32] "TNF signaling pathway"
## [33] "JAK-STAT signaling pathway"
## [34] "Virion - Herpesvirus"
## [35] "Fc epsilon RI signaling pathway"
## [36] "Viral protein interaction with cytokine and cytokine receptor"
## [37] "PD-L1 expression and PD-1 checkpoint pathway in cancer"
## [38] "Antigen processing and presentation"KEGG Pathways - Geni DOWN REGOLATI
downgene_entrez <- down_reg$ENTREZID
Kegg_res <- gost(downgene_entrez, sources="KEGG", organism = "mmusculus", significant = TRUE,user_threshold = 0.05)
Kegg_res$result$term_name## [1] "Focal adhesion"
## [2] "HIF-1 signaling pathway"
## [3] "MAPK signaling pathway"
## [4] "Protein processing in endoplasmic reticulum"
## [5] "Fluid shear stress and atherosclerosis"
## [6] "PI3K-Akt signaling pathway"
## [7] "mTOR signaling pathway"
## [8] "Proteoglycans in cancer"
## [9] "Lysosome"
## [10] "Central carbon metabolism in cancer"
## [11] "Metabolic pathways"
## [12] "Choline metabolism in cancer"
## [13] "ECM-receptor interaction"
## [14] "Sphingolipid signaling pathway"
## [15] "Rap1 signaling pathway"
## [16] "Autophagy - animal"
## [17] "Efferocytosis"
## [18] "Thyroid hormone signaling pathway"
## [19] "Renal cell carcinoma"
## [20] "Amoebiasis"KEGG Pathways - Considerazioni finali
“La sovraespressione di pathways immunitari e infiammatori suggerisce una maggiore attivazione della risposta immunitaria, che può essere positiva se il sistema immunitario riesce a riconoscere e attaccare le cellule tumorali. Tuttavia, l’infiammazione cronica può anche favorire la progressione tumorale.
La sottoespressione di pathways legati alla crescita, adesione e metabolismo cellulare suggerisce una ridotta capacità del tumore di proliferare, sopravvivere e metastatizzare, il che è generalmente positivo per il controllo del tumore. Conclusione:
Nel complesso, il silenziamento di BRCA2 sembra avere un effetto misto, ma prevalentemente positivo rispetto al controllo del tumore al seno. La riduzione della proliferazione e della capacità metastatica delle cellule tumorali, insieme a una potenziale maggiore suscettibilità alle terapie immunitarie, suggerisce che il silenziamento di BRCA2 potrebbe limitare la progressione del tumore al seno. Tuttavia, è essenziale considerare il contesto specifico e ulteriori studi sperimentali per confermare questa valutazione.”
Un’ulteriore informazione interessante è che separando up e down regolati otteniamo un numero minore di pathway rispetto al considerarli entrambi contemporaneamente
#install.packages("remotes")
#remotes::install_github("jokergoo/simplifyEnrichment")
#per poter usare la funzione che fa un subset di un enrichResult: subset_enrichResult(elemento, quantità)Step 8: Network Analysis
Similarità dei campioni
Come ultima parte di questo progetto, ci si pone l’obbiettivo di analizzare i nostri campioni in termini di rete per vedere come questi di distribuiscono in termini di similarità, prima di tutto andiamo calcolare la matrice delle distanze della trasposta dei nostri geni
library(igraph)
library(philentropy)
Gene_Expression <- assay(sef, "tmm")
dist <- as.matrix(distance(t(Gene_Expression), method = "euclidean", use.row.names = TRUE))ottendo le seguenti distanze
dist## BRCA2-1 BRCA2-2 BRCA2-3 BRCA2-4 BRCA2-5 WT-1 WT-2
## BRCA2-1 0.000 3950.972 4078.287 5861.096 3395.781 16882.499 11319.844
## BRCA2-2 3950.972 0.000 3082.494 7578.420 4259.664 20038.013 13630.750
## BRCA2-3 4078.287 3082.494 0.000 7603.953 4353.241 19291.317 13088.773
## BRCA2-4 5861.096 7578.420 7603.953 0.000 4801.625 15145.999 9625.304
## BRCA2-5 3395.781 4259.664 4353.241 4801.625 0.000 17161.324 10957.973
## WT-1 16882.499 20038.013 19291.317 15145.999 17161.324 0.000 9676.455
## WT-2 11319.844 13630.750 13088.773 9625.304 10957.973 9676.455 0.000
## WT-3 15141.359 17802.041 17175.606 13476.008 15245.875 6690.066 5241.374
## WT-4 20013.154 22873.268 22292.685 17756.943 20085.051 6072.720 10596.133
## WT-5 20258.509 23139.703 22624.269 17262.267 20062.953 6839.825 11533.335
## WT-3 WT-4 WT-5
## BRCA2-1 15141.359 20013.154 20258.509
## BRCA2-2 17802.041 22873.268 23139.703
## BRCA2-3 17175.606 22292.685 22624.269
## BRCA2-4 13476.008 17756.943 17262.267
## BRCA2-5 15245.875 20085.051 20062.953
## WT-1 6690.066 6072.720 6839.825
## WT-2 5241.374 10596.133 11533.335
## WT-3 0.000 6295.002 8437.072
## WT-4 6295.002 0.000 4684.591
## WT-5 8437.072 4684.591 0.000Rispetto alle distanze ottenute, andiamo a “tagliare” i collegamenti per quelli che cadono al di sotto del quantile calcolato al 50% per non avere una rete densa e per poter identificare i sottogruppi, inoltre la nostra rete non è pesata quindi sostituiamo al valore di distanza 1 per indicare la presenza di un collegamento ed inoltre annulliamo anche la diagonale principale per evitare loop locali
dist[dist<quantile(dist,0.5)] <- 1
dist[dist>=quantile(dist,0.5)] <- 0
diag(dist) <- 0
dist## BRCA2-1 BRCA2-2 BRCA2-3 BRCA2-4 BRCA2-5 WT-1 WT-2 WT-3 WT-4 WT-5
## BRCA2-1 0 1 1 1 1 0 0 0 0 0
## BRCA2-2 1 0 1 1 1 0 0 0 0 0
## BRCA2-3 1 1 0 1 1 0 0 0 0 0
## BRCA2-4 1 1 1 0 1 0 1 0 0 0
## BRCA2-5 1 1 1 1 0 0 0 0 0 0
## WT-1 0 0 0 0 0 0 1 1 1 1
## WT-2 0 0 0 1 0 1 0 1 1 0
## WT-3 0 0 0 0 0 1 1 0 1 1
## WT-4 0 0 0 0 0 1 1 1 0 1
## WT-5 0 0 0 0 0 1 0 1 1 0a questo punto abbiamo una matrice che ci permette di costruire un grafico in maniera leggera in termini computazionali costruendo un oggeto di tipo iGraph
net <- graph_from_adjacency_matrix(dist, mode="undirected")
plot(net, vertex.size=5,vertex.label.cex=0.5, vertex.frame.color="#ffffff",vertex.color="green")
Osservando come si creino due reti ben separate tra i campioni silenziati e quelli non silenziati dando validità in termini di robustezza alle nostre analisi perchè i campioni dei due gruppi sono simili tra loro.
Clustering dell’enrichmentResult del GO - Ovewview
#remotes::install_github("jmw86069/jamenrich")library(multienrichjam)
library(jamba)
library(colorjam)
go_net <- (go_sym)[go_sym@result$ONTOLOGY == "BP", ]
go_net<-enrichMapJam(go_net)
plot(go_net, vertex.size=8,vertex.label="", vertex.frame.color="#ffffff",vertex.color="green")
Possiamo osservare come nella nostra rete siano presenti dei gruppi più o meno grandi e vogliamo indentificarli con l’algoritmo Louvaine, previa una rimozione dei nodi isolati
ind <- which((degree(go_net)) == 0)
go_net<- delete_vertices(go_net, ind)cl <- cluster_louvain(go_net)
plot(go_net, vertex.size=5,vertex.label="", vertex.frame.color="#ffffff",vertex.color=cl$membership) 
Abbiamo quindi identificato i nostri clusters, per capire questi clusters se sono up o down regolati e quindi, nel caso in esame del silenziamento di BRCA2 vanno ad inibire o stimolare determinati processi bisogna considerare l’up e down regolation dei geni e non la semplice espressione differenziale come abbiamo appena fatto
Clustering dell’enrichment Result del GO - Geni UP REGOLATI
upgene_sym <- up_reg$SYMBOL
upgo_sym <- enrichGO(upgene_sym, OrgDb = organism,
keyType ="SYMBOL", ont = "ALL")
upgo_net <- (upgo_sym)[upgo_sym@result$ONTOLOGY == "BP", ]
upgo_net<-enrichMapJam(upgo_net)
ind <- which((degree(upgo_net)) == 0)
upgo_net<- delete_vertices(upgo_net, ind)
cl <- cluster_louvain(upgo_net)
plot(upgo_net, vertex.size=5,vertex.label="", vertex.frame.color="#ffffff",vertex.color=cl$membership) 
Identificando 3 cluster ben distinti che possiamo studiare
df <- data.frame(length, index)
for (i in seq_along(cl)) {
df <- rbind(df, data.frame(length = length(unlist(cl[i])), index = i))
}
df_clust <- df[order(df$length, decreasing = TRUE), ]
df_clust## length index
## 2 25 2
## 1 22 1
## 3 2 3for (i in 1:nrow(df_clust)) {
num <- df_clust[i, 2]
cat("Cluster", num, ": è formato da\n")
cat(paste(cl[[num]], collapse = ",\n"), "\n\n")
}## Cluster 2 : è formato da
## leukocyte cell-cell adhesion,
## regulation of leukocyte cell-cell adhesion,
## regulation of T cell activation,
## positive regulation of leukocyte cell-cell adhesion,
## lymphocyte differentiation,
## positive regulation of cell-cell adhesion,
## regulation of immune effector process,
## lymphocyte proliferation,
## positive regulation of cell activation,
## leukocyte proliferation,
## mononuclear cell proliferation,
## positive regulation of T cell activation,
## positive regulation of leukocyte activation,
## T cell differentiation,
## T cell proliferation,
## positive regulation of lymphocyte activation,
## regulation of leukocyte differentiation,
## regulation of hemopoiesis,
## regulation of leukocyte proliferation,
## regulation of lymphocyte proliferation,
## regulation of lymphocyte differentiation,
## regulation of mononuclear cell proliferation,
## leukocyte activation involved in immune response,
## positive regulation of immune effector process,
## cell activation involved in immune response
##
## Cluster 1 : è formato da
## defense response to virus,
## defense response to symbiont,
## regulation of innate immune response,
## immune response-regulating signaling pathway,
## immune response-activating signaling pathway,
## positive regulation of innate immune response,
## positive regulation of response to biotic stimulus,
## activation of innate immune response,
## response to interferon-beta,
## cellular response to interferon-beta,
## cytokine-mediated signaling pathway,
## innate immune response-activating signaling pathway,
## pattern recognition receptor signaling pathway,
## cytosolic pattern recognition receptor signaling pathway,
## interferon-mediated signaling pathway,
## negative regulation of viral process,
## negative regulation of cytokine production,
## type I interferon production,
## regulation of viral life cycle,
## regulation of inflammatory response,
## antiviral innate immune response,
## response to virus
##
## Cluster 3 : è formato da
## peptidyl-tyrosine phosphorylation,
## peptidyl-tyrosine modificationAbbiamo quindi trovato per ora questi 3 cluster di processi biologici che nel caso del silenziamento di BRCA2 sono sovra-stimolati
Clustering dell’enrichment Result del GO- Geni DOWN REGOLATI
downgene_sym <- down_reg$SYMBOL
downgo_sym <- enrichGO(downgene_sym, OrgDb = organism,
keyType ="SYMBOL", ont = "ALL")
downgo_net <- (downgo_sym)[downgo_sym@result$ONTOLOGY == "BP", ]
downgo_net<-enrichMapJam(downgo_net)
ind <- which((degree(downgo_net)) == 0)
downgo_net<- delete_vertices(downgo_net, ind)
cl <- cluster_louvain(downgo_net)
plot(downgo_net, vertex.size=5,vertex.label="", vertex.frame.color="#ffffff",vertex.color=cl$membership) 
Identificando 9 cluster ben distinti che possiamo studiare
df <- data.frame(length, index)
for (i in seq_along(cl)) {
df <- rbind(df, data.frame(length = length(unlist(cl[i])), index = i))
}
df_clust <- df[order(df$length, decreasing = TRUE), ]
df_clust## length index
## 1 11 1
## 5 6 5
## 3 4 3
## 4 4 4
## 6 4 6
## 2 3 2
## 7 3 7
## 8 3 8
## 9 3 9for (i in 1:nrow(df_clust)) {
num <- df_clust[i, 2]
cat("Cluster", num, ": è formato da\n")
cat(paste(cl[[num]], collapse = ",\n"), "\n\n")
}## Cluster 1 : è formato da
## regulation of vasculature development,
## regulation of angiogenesis,
## epithelial cell migration,
## epithelium migration,
## tissue migration,
## mesenchymal cell differentiation,
## regulation of epithelial cell migration,
## mesenchyme development,
## positive regulation of angiogenesis,
## positive regulation of vasculature development,
## ameboidal-type cell migration
##
## Cluster 5 : è formato da
## regulation of actin filament-based process,
## regulation of supramolecular fiber organization,
## actin filament organization,
## regulation of cellular component size,
## regulation of actin cytoskeleton organization,
## regulation of actin filament organization
##
## Cluster 3 : è formato da
## ossification,
## bone mineralization,
## biomineral tissue development,
## bone development
##
## Cluster 4 : è formato da
## chemotaxis,
## taxis,
## leukocyte migration,
## cell chemotaxis
##
## Cluster 6 : è formato da
## response to peptide,
## cellular response to peptide,
## cellular response to peptide hormone stimulus,
## response to peptide hormone
##
## Cluster 2 : è formato da
## cell-substrate adhesion,
## cell-matrix adhesion,
## regulation of cell-substrate adhesion
##
## Cluster 7 : è formato da
## extracellular matrix organization,
## external encapsulating structure organization,
## extracellular structure organization
##
## Cluster 8 : è formato da
## response to endoplasmic reticulum stress,
## response to topologically incorrect protein,
## response to unfolded protein
##
## Cluster 9 : è formato da
## positive regulation of cell projection organization,
## regulation of neurogenesis,
## gliogenesisAbbiamo quindi trovato per ora questi 9 cluster di processi biologici che nel caso del silenziamento di BRCA2 sono inibiti.
Clustering dell’enrichment Result del GO - Considerazioni finali
Alla luce dei risultati possiamo trarre le seguenti conslusioni riguardo all’analisi dell’arricchimento fatta con la gene ontology:
“Nel complesso, i geni up-regolati mostrano un miglioramento della risposta immunitaria che può essere positivo per la sorveglianza e l’eliminazione delle cellule tumorali, ma possono anche promuovere l’infiammazione cronica che può avere effetti pro-tumorali. D’altra parte, i geni down-regolati riducono processi cruciali per la crescita, la sopravvivenza e la metastasi delle cellule tumorali, suggerendo un effetto complessivamente positivo contro il tumore.”
Inoltre:
“La clusterizzazione non è stata inutile; al contrario, ha fornito un quadro più chiaro e dettagliato dei processi biologici coinvolti. Ha permesso di valutare meglio l’impatto del silenziamento di BRCA2 e di identificare pathways specifici che possono essere target terapeutici.
In sintesi, mentre una lista di geni regolati è utile, la clusterizzazione offre un valore aggiunto significativo per l’interpretazione biologica e la valutazione finale.”
Step : Heatmap Analisys
In questo ultimo step, diamo uno sguardo alle heatmap e osserviamo diverse informazioni interessanti
Analisi generale
Guardiamo in primis un po’ di heatmap dei geni differenzialmente espressi in generale, quelli up regolati e down regolati
library(pheatmap)
pheatmap(assay(sef)[rownames(assay(sef,"tmm")) %in% diff_exp$ENTREZID, ], scale = "row", cluster_rows = TRUE, cluster_cols = TRUE,
show_rownames = FALSE, show_colnames = TRUE,
main = "Heatmap dei Geni Differenzialmente Espressi")
pheatmap(assay(sef)[rownames(assay(sef,"tmm")) %in% up_reg$ENTREZID, ], scale = "row", cluster_rows = TRUE, cluster_cols = TRUE,
show_rownames = FALSE, show_colnames = TRUE,
main = "Heatmap dei Geni Up_regolati")
pheatmap(assay(sef)[rownames(assay(sef,"tmm")) %in% down_reg$ENTREZID, ], scale = "row", cluster_rows = TRUE, cluster_cols = TRUE,
show_rownames = FALSE, show_colnames = TRUE,
main = "Heatmap dei Geni Down_regolati")
Osserviamo come ci sia una netta separazione in termini di gruppi tra BRCA e WT ad eccezione della mappa dei down-regolati in cui il campione WT-4 viene clusterizzato nel gruppo dei BRCA mentre il campione BRCA2-5 viene clusterizzato nel gruppo dei WT
Heatmap in un pathway
Dato che si parla di cancro, uno dei pathway sicuramente interessante è “immune response-regulating signaling pathway” identificato nei GO-Terms che risulta up-regolato
upgene_sym <- up_reg$ENTREZID
upgo_sym <- enrichGO(upgene_sym, OrgDb = organism,
keyType ="ENTREZID", ont = "BP")
pos<-which(upgo_sym@result$Description == "immune response-regulating signaling pathway")
pathway_genes <- strsplit(upgo_sym@result[pos,]$geneID, split = "/")[[1]]
pheatmap(assay(sef, "tmm")[rownames(assay(sef)) %in% pathway_genes, ], scale = "row", cluster_rows = TRUE, cluster_cols = TRUE,
show_rownames = FALSE, show_colnames = TRUE,
main = "Heatmap dei Geni Up Regolati del Pathway")
Possiamo osservare come, in questo pathway, tutti i geni contribuiscano allo stesso modo perchè non ci sono regioni orizzontali particolamente intense rispetto alle altre (nella partizione BRCA2) anzi, le regioni di maggior intensità per ogni campione di distribuiscono su sotto-gruppi di geni diversi, questo suggerisce che il Pathway di interesse si verifichi effettivamente piuttosto che, come nel caso del pathway del covid 19 questo venga identificato in maniera erronea